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ABSTRACT 

We present the first three-dimensional fully kinetic electromagnetic relativistic particle-in-cell simula- 
tions of the collision of two interpenetrating plasma shells. The highly accurate plasma-kinetic "particle- 
in-cell" (with the total of 10 8 particles) parallel code OSIRIS has been used. Our simulations show: (i) 
the generation of long-lived near-equipartition (electro)magnetic fields, (ii) non-thermal particle accel- 
eration, and (iii) short-scale to long-scale magnetic field evolution, in the collision region. Our results 
provide new insights into the magnetic field generation and particle acceleration in relativistic and sub- 
relativistic colliding streams of particles, which are present in gamma-ray bursters, supernova remnants, 
relativistic jets, pulsar winds, etc.. 

Subject headings: magnetic fields - instabilities - acceleration of particles 



sic scenario unstable to electromagnetic and/or electro- 
static plasma instabilities. To probe the full nonlinear 
dynamics and the saturated state of this system it is nec- 
essary to employ kinetic numerical simulations. We use 
the fully electromagnetic relativistic particle-in-cell (PIC) 
code OSIRIS (Hemker, 2000; Fonseca et al, 2002) to per- 
form the first three-dimensional kinetic simulations of the 
collision of two plasma shells, and to observe the three- 
dimensional features of the electromagnetic filamentation 
instability, or Weibel instability. We compare the collision 
of weakly-relativistic and ultra-relativistic plasma shells. 
We also point out the most distinct features of the three- 
dimensional simulations, not present in lower dimensional 
simulations. Finally, we discuss the key points raised by 
this work and the directions for future work. 

2. SIMULATION MODEL 

In this paper we illustrate the main features of the col- 
lision of two plasma shells, using the PIC code OSIRIS 
(Hemker, 2000; Fonseca et al., 2002). A general description 
of PIC codes is presented by Dawson (1980) and Birdsall 
& Langdon (1982). This scenario is pervasive in astro- 
physical scenarios, and it is of particular relevance in the 
early formation stages of collisionless shocks. Our simula- 
tions are directly relevant to internal shocks in gamma ray 
bursts, in connection with collisions of electron-positron 
fireball shells. 

The simulations were performed on a 256 x 256 x 100 
grid (the axes are labeled as xl, x2, x3) with a total of 105 
million particles for 2900 time steps, with periodic bound- 
ary conditions. A parameter scan was done in order to 
guarantee the grid resolution and the size of the simula- 
tion box do not affect the simulation results. Temporal 
and spatial scales in the simulations are normalized to the 
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1. INTRODUCTION 

Various violent astrophysical phenomena, such as 
gamma-ray bursts (GRBs), supernova explosions, pulsar 
wind outflows, knots in relativistic jets, are known (or 
believed) to be rather strong sources of synchrotron ra- 
diation and cosmic rays, indicating the presence of near- 
equipartition magnetic fields and strong particle acceler- 
ation. These conditions are, in general, associated with 
shocks or regions with colliding streams of particles. How 
the magnetic fields are generated, how long do they live 
and how particles are accelerated arc still open questions, 
which can only be definitely addressed via fully kinetic 
three-dimensional numerical simulations. Most astrophys- 
ical shocks are collisionless because dissipation is domi- 
nated by particle scattering in turbulent electromagnetic 
fields rather than particle-particle collisions (Sagdeev, 
1966). Plasma instabilities driven by streaming particles 
are responsible for the excitation of these turbulent electro- 
magnetic fields. The Weibel instability (Weibel, 1959) has 
received considerable interest as a possible robust mech- 
anism for the production of sub-equipartition long-lived 
magnetic fields and energetic particles in GRBs (Medvedev 
& Loeb, 1999; Pruet et al., 2001) and pulsar winds (Kaz- 
imura et al., 1998). However, Gruzinov (2001a) put the 
simulation results of Kazimura et al. (1998) into question, 
claiming that the produced fields cannot survive far behind 
the shock. Previously, Arons and co-workers also pointed 
out the role of the Weibel instability in the downstream 
region of relativistic magnetosonic shocks (Langdon et al., 
1988; Gallant et al., 1989). 

In this Letter, we examine the dynamics of two col- 
liding inter-penetrating unmagnetized plasma shells with 
zero net charge. This is the most simple model for the 
collision region of two plasma shells, as well as a clas- 



2 



inverse electron plasma frequency to p } = (47re 2 n e /m e ) 1//2 , 
and the collisionless skin depth A e = c/uj pe , where n e is 
the electron-positron shell density; charges and masses are 
measured in the units of the electron charge e and the elec- 
tron mass m e . In these normalized units, the box size is 
V = 25.6 x 25.6 x 10.0 (c/ui pe ) 3 , and the simulations ran for 
150.0 u>~J-. In all runs, energy is conserved down to 0.025 
%. In our simulations, a charge-neutral plasma shell, con- 
sisting of electrons and positrons and moving with a bulk 
momentum u 3 = 70W3/C along the x3 direction (vertical di- 
rection), penetrates into a similar plasma shell, moving in 
the opposite direction with the same momentum u 3 (here 
70 is the initial Lorentz factor of the particles and v 3 is the 
bulk velocity of a shell along x3). Thus, at t = there are 
two groups of particles moving in opposite directions in the 
center of mass frame and occupying the entire simulation 
volume. The particles in both groups (shells) have a ther- 
mal spread with the rms momentum u t h = Jo v th — 0.1. 
The system has no net charge and no net current, and 
initially the electric and magnetic fields are set to zero. 
We performed both sub-relativistic (113 = 0.6, 70 ~ 1.17, 
Vth — 0.085c) and ultra-relativistic (w.3 = 10.0, 70 ~ 10.05, 
Vth — 0.01c) simulations. By using the mass ratio of 1 
(electron-positron plasma), we can follow more clearly the 
crucial features of the instability process. Note that in the 
dimensionless units used, the results of the simulations 
are equally applicable to the collision of proton - anti- 
proton shells as well. We have also performed collisions 
of electron-proton plasma shells and we have observed a 
two-stage instability process, with each stage evolving on 
the electron and proton collisionless time scales, following 
the pattern observed for the electron-positron shells (see 
Medvedcv & Loeb (1999), and Tonge (2002)). 

3. NUMERICAL RESULTS AND PHYSICAL PICTURE 

The fundamental issue to be addressed here is the level 
of the electromagnetic field generated via plasma insta- 
bilities during the collision of two plasma shells, as well 
as the saturated state of the particles and fields. As 
we show, the growth rate is so short that the saturated 
state is all that matters for astrophysical conditions. In 
Figure 1, we present the temporal evolution of the to- 
tal energy in the produced magnetic and electric fields 
normalized to the initial total kinetic energy in the sys- 
tem (e p ), cb = J B 2 dV/8Tre p and e E — J E 2 dV/8ne pi 
for the sub-relativistic (1*3 = 0.6) and ultra-relativistic 
(113 — 10) scenarios. In our normalized units, the ini- 
tial total kinetic energy of particles in the simulations 
is e p = 4 x (70 - 1) x 10.0 x (25.6) 2 . During the lin- 
ear stage of the instability we observed the rapid gen- 
eration of a strong magnetic field, which predominantly 
lies in the xla;2-plane, i.e., perpendicular to the direction 
of motion of the plasma shells. The magnetic field en- 
ergy density reaches ~ 5% of e p for u 3 = 0.6 and ~ 20% 
for 113 = 10.0. In all cases, the produced electric field 
is significantly weaker than the magnetic field. After the 
linear stage, the instability saturates and the energy in 
the magnetic field decays rapidly (on the collisionless time 
scale w~ e ), until it reaches a quasi-steady level with no or 
very slow decay on a time scale much longer than uj p< }. 
These features are present in both the weakly relativistic 
and ultra-relativistic conditions. Such long-lived magnetic 



fields, associated with the Weibel instability, have also 
been observed in laboratory experiments (Lai et al., 1997). 
The linear growth rate agrees well with the theoretical es- 
timates for the full electromagnetic instability (maximum 
growth rate, T, given by {T/uj pe ) 2 ~ 2/3 2 /7 (l + f3 t h) 2 , 
for 70 > 1, and /3 th < 1) (Silva et al., 2002). Figure 
1 also illustrates the main features of the Weibel insta- 
bility, in particular, the linear growth rate scaling with 

— 1/2 

7 and the scaling of the saturation level of the mag- 
netic field with u^ 1 . The saturation level of the magnetic 
field, £> sat , i.e., when B-field growth stalls, can be de- 
termined by combining the analytical results for T, and 
an estimate for the bounce frequency of the particles, 
^bounce, in the magnetic wells generated by the Weibel 
instability. Saturation occurs when T ~ ^bounce ? 

yield- 

I 1/2 1/2 1/2 

ing B sa t ~ V8irm p ' nj (3 c-f /(l + Ah), where m p (q p ) 
is the mass(absolute charge) of the particles in the cloud. 
Electron-proton clouds will drive higher levels of saturated 
B-field, by a factor of (m p /m e ) 1/2 (Tonge, 2002), but on 

a longer time scale since T oc m p 1 ^ 2 . However, cb is inde- 
pendent of n e , 70 and m p {q p ) (see also Medvedev & Loeb 
(1999) for a discussion). Thermal effects are known to play 
an important role in the evolution of the Weibel instability 
(Silva et al., 2002). Higher plasma temperatures will lead 
to a decrease of T and B sat , and longer transverse wave- 
length filaments, within the orders of magnitude observed 
in our simulations (Silva et al. in preparation). The final 
steady-state level of the magnetic field is still quite high, 
up to 0.25% of the initial total kinetic energy in the sub- 
relativistic case and an order of magnitude higher in the 
ultra-relativistic conditions. This residual magnetic field 
is mostly perpendicular, i.e., the xl and x2 components 
are dominant. 

Physically, the inhomogeneities in the current will gen- 
erate inhomogeneities in the magnetic field, which in turn 
will enhance the current inhomogeneities (thus closing the 
instability feedback loop), and generating a large num- 
ber of current filaments with oppositely directed currents. 
The magnetic field associated with these currents has also 
a filamentary structure, as illustrated in Figure 2. In the 
early stages of the instability, one observes randomly dis- 
tributed current and magnetic filaments. Figure 2a shows 
small-scale tilted iso-levels of the magnetic field energy 
density. When the instability enters the saturation phase 
(t w 15 — 30 u p e), current filaments begin to interact 
with each other, forcing like currents to approach each 
other and merge. During this phase, initially randomly ori- 
ented filaments cross each other to form a more organized, 
large-scale quasi-regular pattern. The strong decrease in 
the magnetic field energy is associated with a topological 
change in the structure of currents and magnetic fields, cf. 
Figure 2a, 2b, and 2c. After saturation (t ;> 30 u>~^), the 
filament coalescence continues, as indicated by the increase 
of the correlation scale, ~ fc _1 , of the .B-field in Figure 3. 
However, the spatial distribution of currents is now quite 
regular, so that filaments with opposite polarity no longer 
cross each other but simply interchange, staying always far 
away. The total magnetic field energy does not change any 
more. Note that the residual magnetic field is highly inho- 
mogeneous, seen as a collection of magnetic field domains 
or " bubbles" . The amplitude of the field in the bubbles 
is close to equipartition. Therefore, the overall decrease 
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of the S-field energy is mostly associated with decreas- 
ing filling factor of the field. Note also that the magnetic 
domains separate current filaments of opposite polarity. 
The discussed temporal evolution can be easily followed in 
the magnetic field energy spectral density, Figure 3. The 
small-scale magnetic field is generated in the linear stage 
of the instability, attains its maximum value, and after sat- 
uration it evolves to larger length scales (i.e., smaller |fc|) 
with the characteristic e- folding time ~ ub = eB/m e c 
(Medvedev, in preparation); note that us/^pe = V^b 
and this is valid for relativistic particles as well. 

The topological evolution of the magnetic field is ac- 
companied by strong heating and non-thermal particle ac- 
celeration, as illustrated in Figure 4. Generation of non- 
thermal fast particles is more pronounced in the highest 
£>-field scenarios. We have observed thermal (rms) mo- 
menta u t h as high as the initial bulk momentum u 3 . The 
magnetic field energy grows in the early stage by slow- 
ing down the plasma shells. Saturation is achieved by 
the combination of transverse energy spread and near- 
cquipartition magnetic field generation. After saturation, 
the energy stored in the magnetic field is transferred back 
to the plasma particles, leading to strong heating and the 
generation of high energy tails in the distribution function, 
with energies up to 4(70 — 1) for both sub-relativistic and 
ultra-relativistic plasma shells. The presence of such high- 
energy particles is fundamental to provide the mildly rela- 
tivistic particles to be injected in other accelerating struc- 
tures, such as collisionless shocks. Comparison of Figures 
3 and 4 also shows that non-thermal particle acceleration 
and strong plasma heating are correlated with the topo- 
logical evolution of the magnetic field when it evolves from 
small-scale structures to large-scale structures. 

To elucidate the difference between the evolution in 
three dimensions vs. two dimensions, we performed several 
2D runs using the same code: the runs in a plane paral- 
lel to the streaming velocity lead to an under-estimation 
of the magnetic field energy, while 2D runs in the plane 
perpendicular to the streaming direction of motion show 
suppression of the excitation of relativistic plasma waves 
along this direction and, hence, overestimate the magnetic 
field energy. 

A critical question for the validity of the GRB syn- 
chrotron shock model (see, e.g., a review by Piran, 1999) 
is whether the magnetic field produced in the collision re- 
gion survives for a sufficiently long time (e.g., compara- 
ble to the synchrotron cooling time), typically ~ 10 w~ e 
for prompt GRB emission and early afterglow and up to 
<~ lO 10 ^,, 1 for late (radio) afterglow, as measured in the 
collision frame. Our simulations show very weak or no 
field dissipation on time-scales ~ 150 w~ e in the collision 
region. In principle, on much longer time-scales, the B- 
field may be destroyed via reconnection. The back-of-the- 
envelope estimate of the reconnection time yields: t rec ~ 

{Vreck)- 1 ~ (O.lVAk)' 1 ~ (O.le^fcAe)" 1 ~ 100(fcA e )-\ 



where v rec ~ O.Iva is the reconnection velocity (Biskamp 
& Schwarz, 2001) and va = B/ (47TTO e n e ) 1 / 2 is the electron 
Alfven speed. Since k decreases with a shorter e-folding 
time ~ 1/y/es ~ 10, reconnection slows down at large 
times and the field is able to survive for quite a while. 
In addition, a fundamental three-dimensional feature is 
the small tilting of current filaments. This is intrinsically 
three-dimensional and it is a manifestation of the coupling 
between the electromagnetic modes (k in the xlx2-plane) 
and the electrostatic modes (k parallel to the streaming 
velocity of the clouds), making magnetic field reconnec- 
tion much more complicated (Silva et al. in preparation). 

4. CONCLUSIONS 

We presented the first self-consistent three-dimensional 
simulations of the fields present in collisions of plasma 
shells where the electromagnetic two-stream (or Weibel) 
instability develops. Our results demonstrate that this 
instability in three-dimensions is able to generate sub- 
equipartition quasi-static long-lived magnetic fields on the 
collisionless temporal and spatial scales in the collision re- 
gion, giving credence to the predictions of Medvedev & 
Loeb (1999) and Pruct ct al. (2001). After the linear stage 
of the instability, we first observe the decay of the mag- 
netic field energy, as also observed by Gruzinov (2001b), 
followed by the evolution to a residual saturated magnetic 
field energy density. These fields maintain a strong satu- 
rated level on time-scales much longer than collisionless, 
at least for the duration of the simulations. The next 
step is to increase the simulation region spatial dimen- 
sion in order to determine the spatial spread of the gen- 
erated field. The obtained values of the equipartition pa- 
rameter (cb ~ 3 x 10~ 2 — 3 x 10~ 3 for ultra- and sub- 
relativistic shocks) agree well with the values of in- 
ferred from GRB afterglows (Panaitescu & Kumar, 2002). 
Furthermore, we show that strong transverse temperature 
increase and non-thermal particle acceleration occur when 
the instability saturates. The generated magnetic field 
evolves from small wavelengths to long wavelengths. Such 
a behavior may explain the observed evolution of the soft 
spectral index a in time-resolved GRB spectra (Ghirlanda 
et al., 2002; Medvedev (in preparation)) as a jitter-to- 
synchrotron spectral transition (Medvedev, 2000). Our 
results indicate that the fields necessary in the early for- 
mation stages of a shock front can be easily created via 
plasma instabilities of streaming plasmas. These simu- 
lations open the way to the full three-dimensional PIC 
modeling of relativistic collisionless shocks, which neces- 
sarily involve a different simulation geometry, and will be 
presented elsewhere. 
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Fig. 1. — Temporal evolution of the total energy in B± = (B 2 + B 2 ) 1 ' 2 and in E%, normalized to t p . Energy on other components, i.e., in 
B-i and E ± = (E 2 + E 2 ) 1 / 2 is much smaller. 
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Fig. 2. — Current filamcntation: magnetic field energy density in weakly relativistic scenario (70 = 1.17): (a) slightly before saturation, 
t = 13.52 u>pe , (b) after saturation t = 41.60 u>pe , for an iso-surface at 20 % of the maximum energy density; in (c) mass density after 
saturation t = 41.60 ujp e for an iso-surfacc at 75 % of the maximum electron density for species with initial positive current along x3 (73 > 0) 
(red iso-surface) and j'3 < (blue iso-surfacc). Color bar in (a) and (b) represents magnetic field energy density (in simulation units) - plotted 
only values above the isosurface level (@ 20% of maximum). 



Fig. 3. — Temporal evolution of the magnetic field energy spectral density distribution, normalized to the peak spectral density (70 = 1.17). 



Fig. 4. — Temporal evolution of the energy distribution of plasma particles (color scale = absolute value of charge density in the simulation 
units). 
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